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Resumen 


En esta contribución se abordan diferentes modelos agregados, descritos 
en la literatura especializada para el estudio de descargas de acuíferos, 
incluyendo un modelo de tipo no lineal propuesto por el autor. Con el 
objetivo de estimar los parámetros involucrados, cada modelo presentado 
fue acoplado con el algoritmo metaheurístico de búsqueda global Cuckoo 
Search (CS), el cual genera la aproximación inicial del algoritmo de 
búsqueda local FMINCON de MatLab. Con base en datos reales medidos en 
el acuífero Ciego-Morón, los parámetros que caracterizan el agotamiento 
del acuífero se estimaron de forma satisfactoria, así como el volumen de 
recarga ocasionada por la precipitación. Los resultados obtenidos se 
sustentaron por las elevadas correlaciones entre simulaciones y 
observaciones, demostrando la factibilidad del abordaje expuesto en el 
presente artículo. 


Palabras clave: descargas, acuíferos, estimación de parámetros. 


Abstract 


In this contribution have been approached different lumped models, 
described in the specialized literature for the study of aquifer discharges, 
including a non-linear model proposed by the author. In order to estimate 
the parameters involved, each model presented has been coupled with the 
global search metaheuristic algorithm Cuckoo Search (CS), which 
generates the initial guess for the local search algorithm FMINCON of 
MatLab. Based on real measured data in the aquifer Ciego-Morón were 
successfully estimated the parameters that characterize aquifer's decay and 
recharge volume due to precipitation. The results were supported by the 
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high correlation between simulations and observations, demonstrating the 
feasibility of the approach considered here. 


Keywords: Discharge, aquifers, parameters estimation. 
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Introducción 


El escurrimiento proveniente de corrientes superficiales y manantiales —en 
ausencia de precipitaciones— es alimentado por las reservas subterráneas 
de los acuíferos, las cuales se han acumulado en dichas estructuras 
geológicas en el transcurso de una recarga determinada (Pérez-Franco, 
1995). 


Este proceso puede ser simulado a partir de modelos matemáticos 
(agregados o distribuidos), que establecen relaciones entre la precipitación 
y el escurrimiento. Este tipo de modelos se ha venido utilizando para 
evaluar los recursos hídricos renovables de un sistema subterráneo, 
dependiendo de la complejidad y finalidad de los modelos de gestión 
adoptados. Uno de los elementos que soporta la elección de los modelos de 
referencia es la disponibilidad de datos. No siempre se cuenta con datos 
piezométricos y de explotación que permitan utilizar modelos distribuidos, 
por lo que en varias ocasiones, a pesar de la variabilidad hidráulica del 
acuífero en estudio, debe recurrirse razonablemente a modelos agregados, 
los cuales permiten, utilizando pocos parámetros, globalizar el 
comportamiento del acuífero a partir de datos en ríos o manantiales 
relacionados con las recargas y/o extracciones con cierta precisión. 


De acuerdo con Arenillas et al. (s.f.), estos modelos agregados de 
precipitación-escurrimiento (MPE) forman parte de un grupo específico de 
herramientas computacionales para la simulación hidrológica, cuyo principal 
objetivo es la generación (o reconstrucción) de series temporales de 
escurrimiento, partiendo de una serie de datos de precipitaciones y 
extracciones por bombeo, así como temperatura existente en la zona de 
estudio. 


En dichos modelos intervienen parámetros físicos que son característicos 
del sistema que se pretende simular y, en sentido general, son 
desconocidos, resultando de suma importancia abordar la formulación de 
un problema inverso que conduzca a una estimación óptima de los mismos, 
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pues intervienen en la evaluación de los recursos hídricos del sistema en 
cuestión. 


La solución de un problema inverso puede ser hecha a través de algoritmos 
de optimización, con el propósito de minimizar (o maximizar) cierta función 
objetivo. Básicamente existen dos clases de métodos de optimización: 
determinísticos y estocásticos. Los primeros dependen del gradiente de la 
función objetivo, haciendo uso de las particularidades del funcional; en 
cambio, los métodos estocásticos trabajan con la selección de 
configuraciones aleatorias; por ejemplo, técnicas computacionales 
inspiradas en la naturaleza, con el objetivo de barrer (estocásticamente) y 
el espacio de búsqueda completo, al procurar el óptimo global (Silva-Neto 
8: Becceneri, 2009). 


En el presente artículo se abordarán diferentes modelos agregados para el 
estudio de descargas de acuíferos. Para establecer una comparación entre 
ellos, el autor inicialmente unificará el enfoque teórico que permitirá 
obtener las respectivas expresiones que estimarán los volúmenes de 
descarga en acuíferos kársticos. En la sección 2, al inicio de este estudio, 
será descrito el modelo conceptual MEDA, propuesto por Iglesias (1984), el 
cual aplica la fórmula propuesta por Maillet (citado por Pérez Franco, 1995) 
para el estudio del agotamiento del escurrimiento proveniente de un 
acuífero. Otras formulaciones, como las desarrolladas por Tisson et al., así 
como Forkasiewicz y Paloc (citados por Pérez-Franco, 1995), también serán 
objeto de análisis, con el propósito de formular modelos que permitan 
estimar las descargas en cuestión. En la sección 3, el autor presentará y 
desarrollará un modelo no lineal. La sección 4 describe el problema inverso 
donde es formulada la función objetivo y los algoritmos de búsqueda global 
utilizados, mientras que en la sección 5 se describe el caso de estudio. 
Posteriormente se muestran los resultados obtenidos y su discusión. Por 
último se tienen las conclusiones derivadas de los objetivos perseguidos en 
esta contribución. 


Modelos agregados para cuantificar caudales de descarga 
de acuíferos 


El caudal de determinada fuente, o la descarga de un acuífero en sentido 
general, no se correlaciona con la precipitación con carácter mensual. Es 
decir, las descargas dependen del volumen almacenado en el acuífero y 
éste guarda correlación con las precipitaciones precedentes en varios 
meses e incluso en años anteriores. En cambio, si se consideran la 
precipitación de cierto periodo y el volumen de descarga acontecido a lo 
largo del tiempo como consecuencia de esa precipitación es posible 
establecer una clara interrelación. El método para determinar el volumen 
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de descarga consecuente a una precipitación dada ha sido abordado en 
detalle por Castillo (1981), citado por Arenillas et al. (s.f.). Con base en 
series de precipitación (P) y volúmenes observados (V), se han obtenido 
ecuaciones de ajuste de tipo lineal y potencial, con coeficientes de 
correlación (r) del orden de 0.9 en la forma: 


e Ecuación lineal: 

V=mP+n (1a) 
e Ecuación potencial: 

V=mP” (1b) 


siendo m y n = parámetros de ajuste. En función del tipo de ecuación 
adoptada, los parámetros m y n presentan dimensiones. Esta observación 
es pertinente, pues en Arenillas et al. (s.f.), ambos parámetros son 
literalmente presentados como adimensionales, lo cual constituye, sin 
duda, un error en dicha publicación. En efecto, para el enfoque lineal (m) = 
L? y (n) = L?, mientras que para el potencial (m) = L?* ” y (n) = 
adimensional, al tratarse del exponente. 


Formulación de Iglesias 


Tal y como fue enunciado con anterioridad, Iglesias (1984) introduce el 
modelo propuesto por Maillet (citado por Pérez Franco, 1995) para la curva 
de agotamiento de la descarga de un acuífero, el cual se expresa como: 


Q, = Quer“ (2) 


siendo Q: = descarga del acuífero en el instante de tiempo t; Q. = descarga 
al inicio del periodo analizado; a = coeficiente de agotamiento, 
característico del acuífero. Acorde con Pérez-Franco (1995), este parámetro 
a = 1/7, siendo r el tiempo necesario para que Q; = Qo/e = 0.3678 Qo. 


Tal y como se muestra en la Figura 1, en el ¡-ésimo mes (mes /), el caudal 
de descarga se denota por Q;. En caso de no existir recarga como 
consecuencia de precipitaciones, entonces el caudal del mes ¡ + 1 será 
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denotado como Q;,,. Si se reporta la existencia de una precipitación P;, el 
caudal de descarga sería obviamente superior y se identificará como Q;,;. 


Función objetivo: f(x),x = (X,,X2, «.-, Xa); 

Generar una población inicial de n nidos; 

While (t<MaxGeneration) or (criterio de parada) 

Obtener un cuco aleatoriamente (por ejemplo, /) y se reemplaza esta 
solución introduciendo los vuelos de Lévy; 

Evaluar su calidad/ajuste F; 

Escoger un nido entre n (por ejemplo, 7) aleatoriamente; 

if (E, > F) 

Reemplazar j por la nueva solución; 

end if 

Una fracción (pa) de los peores nidos es abandonada y se construyen 
nuevos; 

Mantener las mejores soluciones/nidos; 

Clasificar las soluciones/nidos y encontrar las mejores; 

Pasar las mejores soluciones a la próxima generación. 

end while 


Figura 1. Pseudocódigo del algoritmo CS. Adaptado de Yang y Deb (2009). 


La precipitación de referencia P;,, acorde con Iglesias (1984), habrá de 
incrementar el volumen almacenado en el acuífero en una cantidad 
equivalente al área comprendida entre las curvas de agotamiento 
representadas en la Figura 2 (Arenillas et al., s.f.). 


Precipitación 


Pi 


Curva de agotamiento para 
Qi+1 


Curva de agotamien 
para Qi 


mes ¡ mes ¡+1 


y= f omar- eva 


mes i+1 mes i+1 
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Figura 2. Curvas de agotamiento en ausencia y presencia de precipitación. 


En efecto, aplicando la ecuación (2), el área bajo las curvas de agotamiento 
identificada como el volumen de descarga V puede ser expresada por la 
diferencia de las siguientes integrales: 


En ES id A Qísa _ Qi qn 
V=V,-V,=f, Qirre “dt f, Qi e dt = 42-44 = 2-4 (8) 
siendo At; = intervalo entre el mes ¡ + 1 y mes /. Como se pretende 


determinar Q;+1, cualesquiera de las expresiones (1) puede ser utilizada. 
Arenillas et al. (s.f.) recomiendan el empleo de la ecuación (1b) 
correspondiente a la estructura potencial de forma generalizada, ya que la 
ecuación (1a) tiene un rango de validez restringido, y puede llegar a dar 
resultados anómalos e incluso negativos. Para evitar dicho inconveniente y 
dar continuidad al desarrollo teórico aquí presentado, el autor será 
consecuente con dicha recomendación, sin ser excluyente con el posible 
empleo de la ecuación (1a). 


Sustituyendo la ecuación (1b) en (3) y despejando Q;+1 es posible obtener 
la siguiente expresión: 


Qisa = (mPP) « +Qje et (4) 


La ecuación (4) es la base conceptual del Modelo para el Estudio de 
Descargas de Acuíferos (MEDA) del Instituto Tecnológico GeoMinero de 
España, y con el mismo es posible realizar las estimaciones pertinentes de 
los caudales de descarga de cada mes, sobre la base del conocimiento de 
precipitación, y caudales del mes precedente y los parámetros m, n y Qu. 
Posteriores análisis del modelo han aportado mejoras notables en su 
formulación conceptual. Entre éstas se puede mencionar la introducción del 
concepto de precipitación útil (P,), definida según Arenillas et al. (s.f.), 
como la diferencia entre la precipitación (P) y la evapotranspiración real 
(Es), la cual forma parte de la recarga en el acuífero y está estrechamente 
relacionada con la temperatura sobre la superficie del suelo (T) y situación 
geográfica, entre otros parámetros del clima. Arenillas et al. (s.f.), al 
analizar valores de evapotranspiración potencial y real en España, 
presentan la siguiente ecuación de cálculo de la precipitación útil para el ¡- 
ésimo mes: 


P,, =P, -— Tf (5) 


donde el parámetro b = exponente, que puede variar entre 1.3 y 1.6. 
Expresiones similares a la anterior pueden ser obtenidas para condiciones 
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locales, características de regiones y climas concretos. No obstante, de 
acuerdo con Arenillas et al. (s.f.), la propuesta abordada no resta 
generalidad al modelo MEDA y, por lo tanto, la ecuación (4) podrá ser 
rescrita como: 


Qis1 = (mP) «+00 (6) 


De igual forma, es posible incluir el volumen de los bombeos extraídos 
mensualmente en el acuífero (W). Esta inclusión, tal y como reportan 
Arenillas et al. (s.f.), presenta la limitación, en ocasiones insalvable, de que 
éstos no produzcan afectaciones de carácter dinámico al punto de 
descarga. Por tanto, de ser considerado el volumen bombeado en el ¡- 
ésimo mes (W;), la ecuación (6) podrá ser replanteada como: 


Oh = (m PR M1) ec +Qyenros 0) 


La ecuación (7), para las consideraciones establecidas hasta su 
conceptualización, representa la expresión más general y utilizable del 
modelo MEDA (de tipo MPE), la cual, en ausencia de datos de temperatura 
y/o bombeos mensuales, presenta como casos particulares las ecuaciones 
(4) y (6), respectivamente. 


Formulación de Tisson, Werner y Dunsquit 


Pérez-Franco (1995) ha desarrollado la siguiente expresión para el 
agotamiento de la descarga de un acuífero: 


Q; == (8) 


- (1+a t)?2 


De forma análoga a la propuesta de Maillet, empleada por Iglesias (1984) 
en el desarrollo de MEDA, el parámetro a se identifica como el coeficiente 
de agotamiento del acuífero. Como ya es conocido « = 1/7; no obstante 
esta vez res el tiempo necesario para que Q: = Q./4. Este razonamiento 
sugiere que para el mismo valor de a, las expresiones (2) y (8) arriban a 
caudales de descarga diferentes, siendo mayores aquellos calculados a 
través de la ecuación (2). Para obtener el volumen de descarga V, con el 
empleo de la ecuación (8) se procede de forma análoga al cálculo de las 
siguientes integrales: 
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o Qi41 co Qipy —Qiga — Qiga — Qiga Q; 
le (1+a yA =S, (1+a t)2 dt = x xa x(1+a At¡)2 (9) 


Tomando en cuenta las mejoras introducidas por Iglesias (1984) (citado 
por Arenillas et al., s.f.) en la cuantificación de V, es posible obtener un 
modelo para la estimación de las descargas de un acuífero (ver ecuación 
(10)), donde es de destacar la alteración que sufre el segundo término de 
la ecuación (7), que caracteriza el agotamiento del acuífero a partir del ¡- 
ésimo mes, como consecuencia de introducir la formulación (8): 


Qis1 = (m Pi —W,) « + (10) 


EEE 
Formulación de Forkasiewicz y Paloc 


Esta fórmula fue propuesta sobre la base de observaciones en una fuente 
kárstica y presenta la siguiente estructura (Pérez-Franco, 1995): 


(11) 


ja Q 
0: = === 
qe 1+6802t 


donde KB = coeficiente de decrecimiento del caudal de descarga. Las 
unidades de este parámetro son L”*-T. Su interpretación puede llevarse a 
cabo igualando la ecuación (11); por ejemplo, con la ecuación (2). Al 
establecer las transformaciones pertinentes es posible encontrar la 
siguiente relación: 


e2at_4 


Q5t ELA 


fB = 


La ecuación anterior no sólo sugiere que KB depende de la capacidad de 
agotamiento del acuífero, expresada a través de «a, sino también del caudal 
de descarga Qo. y además del tiempo. Puede obtenerse un resultado similar 
al emplear las ecuaciones (11) y (8). Sin embargo, su desarrollo no 
reportaría ningún elemento cualitativo adicional al obtenido por la vía aquí 
expuesta. Precisamente, en esta multidependencia radica la complejidad de 
la formulación analizada. Aunque en Pérez-Franco (1995) no se presentan 
elementos adicionales que permitan esclarecer aún más la naturaleza de 
dicho parámetro, es de suponer que sería útil para reconstruir series 
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históricas de descargas de acuíferos kársticos en determinado periodo, y 
por tal motivo ha sido introducida en el presente artículo. El volumen de 
descarga V en este caso podrá expresarse como: 


Y = [7 dt — $? dt = ¿(A ) (13) 


o |] O, 
1+BQÍ,,t 1+B(03,,)+ ds Qi+1 


El modelo que resulta para la estimación de descargas de un acuífero a 
partir de la ecuación (11) presenta una estructura diferente a las 
presentadas en las ecuaciones (7) y (10). Este modelo (ver ecuación (14)) 
aporta además la presencia de un nuevo parámetro (fB), que caracteriza el 
decrecimiento de la descarga de una fuente existente en un acuífero, que 
como ya se ha explicado, es de naturaleza compleja y su determinación 
debe ser rigurosa para no incurrir en errores no deseados que afecten la 
calidad de las estimaciones efectuadas en la descarga. 


Qi (14) 


Qi = Y 7 
[1+807 At;-(m pz -w;)£0; 


Propuesta de modelo no lineal para determinar la descarga 
para un acuífero 


Tal y como se ha podido apreciar en las secciones anteriores, existen varios 
modelos para describir la dinámica del agotamiento de la descarga de un 
acuífero. Uno de los más usados es el propuesto por Maillet (citado por 
Pérez-Franco, 1995). Desde el punto de vista de su estructura (ver 
ecuación (2)), dicho modelo es solución de la siguiente ecuación diferencial 
ordinaria: 


qq (15) 


dt 


sujeta a la siguiente condición inicial Q(0) = Qo. Independiente de los 
niveles de precisión que puede brindar este tipo de formulación, una 
pequeña modificación en la ecuación (15) puede ser introducida, 
considerando que la variación de la descarga en el acuífero es una función 
no lineal. Esta puede ser escrita como: 


IQ (16) 


dt 
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siendo n = exponente y xk = parámetro que caracteriza el agotamiento del 
acuífero, y cuyas unidades son LY T”?. La integración de la ecuación 
(16) conduce a la siguiente solución analítica: 


Q:=[057" — (1 me 6” (17) 
Como ya es usual, el volumen de descarga V podrá ser determinado según: 


A E (18) 
A CS 


Para la actual formulación, resulta de interés para el autor expresar las 
integrales de la ecuación (18) empleando procedimientos apropiados para 
el cálculo de límites, dado que se está en presencia de integrales 
impropias, por lo tanto: 


(1-9) pes E (>) 2- * =- (55) * E 
y tin KM e] E Ra a e 
= =- — lim - 
p6s x(n — 2) k(n-2) too «(y — 2) «(y — 2) 
2-7 * (Q-N) 
— Pla (Qí+1) (19) 


kQ-9) k(Q-9) 


En la ecuación (19), para O < n < 1, los límites del primer y tercer términos 
convergen a cero y el volumen de descarga resulta de la diferencia entre el 
segundo y cuarto término en dicha ecuación. Esta particularidad es 
aprovechada para establecer la formulación del modelo propuesto en esta 
sección, de forma análoga, a como se procedió para obtener las ecuaciones 
(7), (10) y (14). En efecto, el modelo no lineal para establecer las 
descargas de acuíferos, válido para O < n < 1, se expresa por la ecuación 
(20) siguiente: 


dí 


Qirs = [o PE —W,) (2 + [07 — (1 A Ja (20) 


Formulación del problema inverso 
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Acorde con Silva-Neto y Moura-Neto (2005), así como con Silva-Neto y 
Becceneri (2009), los problemas inversos pueden ser clasificados tomando 
en cuenta la dimensión del espacio donde el modelo es definido, así como 
la dimensión del espacio donde los objetos que son estimados pertenecen. 
Esta clasificación se resume en la Tabla 1 y es útil para brindar una idea a 
priori de la dificultad del problema en consideración. 


Tabla 1. Clasificación de problemas inversos. 


Modelo matemático Estimación 
Dimensión finita Dimensión infinita 
(constantes) (función) 
Dimensión finita * Tipo 1 - 
Dimensión infinita ** Tipo 11 Tipo III 


“por ejemplo: sistema de ecuaciones algebraicas; “por ejemplo: sistema de ecuaciones 
diferenciales. 


En los modelos dados por las ecuaciones (7), (10), (14) y (20) es posible 
estimar de modo simultáneo los parámetros que son desconocidos. Como 
la cantidad de mediciones debe ser superior (o al menos igual) a la 
cantidad de incógnitas a ser estimadas (parámetros), su identificación se 
basa en la resolución del problema inverso, tratado como la minimización 
de cierta función objetivo en un espacio de búsqueda de dimensión finita. 


Función objetivo 


La función objetivo usada en este estudio es multicriterio, y ha sido 
adaptada de Barco, Wong y Stenstrom (2008): 


F(8) = 9 (LO 4 y, (PO, y, (O (21) 


max 
Qr 00 7; 


sujeta 0,<0<0,,, siendo O = vector de los parámetros; 8, y 6, son los 
vectores de los respectivos límites inferior y superior de los parámetros; V 
= volumen total de flujo descargado por el acuífero; Q”** = caudal máximo 
de descarga mensual; Q = descarga mensual; el subíndice 'o” denota el /- 
ésimo valor de la serie de caudales observados, mientras que el subíndice '¡ 
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"denota la j¡-ésima simulación del caudal de descarga. Los parámetros w,, 
W> Y W3 representan factores de ponderación (o peso) de cada uno de los 
términos de la función objetivo. La división con respecto a los valores 
observados normaliza dicha función y en la presente contribución será 
asignada igual importancia a cada uno de los términos de la ecuación (21) 
considerando W;, = W2 = wz = 1.0. La cantidad de parámetros escogidos 
para la calibración depende de la estructura de los modelos se presentan 
en las secciones anteriores. De forma general, en la Tabla 2 se presentan 
los posibles parámetros a identificar. 


Tabla 2. Parámetros de los modelos. 


Modelo Vector de parámetros 


Ecuación (7) 
Ecuación (10) 
Ecuación (14) 0=(m,n,b,B)' 
Ecuación (20) 0=(m,n,b,x,n)' 


0=(m,n,b,a)' 


En esta investigación se empleó el algoritmo de búsqueda global Cuckoo 
(Cuckoo search, CS), desarrollado por Yang y Deb (Yang € Deb, 2009; 
Yang € Deb, 2010). Tal y como reportan Civicioglu y Besdok (2011), la 
introducción de este algoritmo puede proveer resultados más robustos que 
algoritmos ampliamente utilizados, como el algoritmo de optimización por 
enjambre de partículas (Particle Swarm Optimization, PSO) desarrollado 
por Clerc y Kennedy (2002), así como el basado en una colonia artificial de 
abejas (Artificial Bee Colony Algorithm, ABC) desarrollado por Karaboga y 
Basturk (2007a, 2007b). Además, se ha optado por emplear un enfoque 
híbrido, acoplando CS con algún método determinístico basado en el 
gradiente. En este caso, será introducida la herramienta FMINCON del 
asistente matemático MatLab, la cual consta de varios algoritmos y 
permite, como bien es conocido por los utilizadores de este asistente, 
encontrar el mínimo de una función no lineal, multidependiente sujeta a 
restricciones. En esta contribución será empleado el algoritmo punto 
interior (interior-point), el cual le incorpora rapidez y precisión al problema 
inverso abordado, al tratarse de minimizar la función objetivo (ecuación 
(21)), sujeta a restricciones de los parámetros. 


El algoritmo CS 
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El algoritmo CS fue desarrollado por Yang y Deb (Yang €: Deb, 2009; Yang 
8 Deb, 2010) basado en la aplicación metaheurística proveniente de la 
naturaleza e inspirado en el comportamiento de pájaros fascinantes, no 
sólo por los sonidos que pueden hacer sino por su agresiva estrategia de 
reproducción. Acorde con dichos investigadores, algunas especies de aves, 
como Ani y Guira cuckoos, colocan sus huevos en nidos comunes, donde 
ellos pueden llegar a remover los huevos de otras especies de pájaros para 
incrementar la probabilidad de producir sus propios huevos. 


Un número considerable de especies de cuckoos (o cucos) presentan un 
comportamiento parásito al momento de anidar, colocando sus huevos en 
el nido de otras especies para que éstas se encarguen de cuidar y alimentar 
a sus polluelos. Pueden llegar a establecerse conflictos directos con los 
cucos intrusos y si una especie descubre que los huevos en su nido no son 
propios, simplemente aparta estos huevos Oo abandona ese nido, y 
construye otro en un lugar diferente. Acorde con Yang y Deb (2009), 
algunas especies de cucos, como el New World brood-parasitic Tapera han 
evolucionado de tal manera que las hembras han podido asemejar en color 
y patrones de los huevos de un número selecto de especies de pájaros. 
Esto reduce la probabilidad de que sus propios huevos sean abandonados e 
incrementa su descendencia. El algoritmo CS emplea las siguientes 
representaciones: 


e Cada huevo en un nido representa una solución y el huevo de cuco 
representa una solución nueva. El propósito es usar las nuevas y 
mejores soluciones (cucos) para reemplazar aquellas que no son tan 
buenas en los nidos. 

e En la forma más simple, cada nido tiene un huevo; sin embargo, el 
algoritmo puede ser extendido a casos más complicados, en los que 
cada nido tiene múltiples huevos, representando un conjunto de 
soluciones. 


El algoritmo CS está basado en tres reglas idealizadas: 


1. Cada cuco coloca un huevo a la vez y lo hace en un nido escogido de 
forma aleatoria. 

2. Los mejores nidos, con alta calidad de huevos, serán transferidos a la 
próxima generación. 

3. El número de enjambres de nidos es predeterminado (Ny) y el huevo 
colocado por el cuco es descubierto por las especies de pájaros con 
una probabilidad pa,= €(0, 1). 


Además, Yang y Deb (2009) descubrieron que el estilo de búsqueda del 
camino aleatorio es mejor ejecutado por los vuelos Lévy. El pseudocódigo 
puede ser resumido como se ilustra en la Figura 2. 


La importante ventaja de este algoritmo es su simplicidad, así como 
también su facilidad para ser implementado computacionalmente. Es válido 
resaltar que el algoritmo propuesto por Yang y Deb (2009) ha sido 
modificado por diferentes investigadores en este campo. El autor del 
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presente artículo ha utilizado el código de CS programado en MatLab por 
Yang, en la Universidad de Cambridge, entre noviembre de 2008 y junio de 
2009, con base en la propuesta original del algoritmo. Como ya se 
mencionó, se ha acoplado a la función FMINCON de MatLab, de manera 
que, en adelante, para hacer referencia a la estrategia adoptada esta 
combinación será denominada CS-FMIN. En la Figura 3 se muestra el 
diagrama de flujo para la estimación de parámetros en modelos de 
descarga de acuíferos. Para la aplicación de CS-FMIN, inicialmente será 
utilizada una versión débil de CS, con el objetivo de generar una 
aproximación inicial para el algoritmo IP. Tal y como muestra la Figura 4, 
en primer lugar, se establece un procedimiento iterativo para estimar los 
parámetros del modelo de descarga analizado, pero con el objetivo de 
mejorar los resultados es que se introduce la estimación definitiva a través 
de FMINCON. 


Aproximación inicial generada por CS 
Estimar vector de parámetros 0 


¿Cumple el criterio de parada? 


STOP 


Figura 3. Diagrama de flujo para la estimación de los parámetros en modelos de 
descargas de acuíferos. 
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Figura 4. Ubicación geográfica y cuenca superficial asociada al área del acuífero Ciego- 
Morón. Adaptado de Vidal-Olivera y González-Abreu (2013). 


Materiales y métodos 


Estudio de caso 


A continuación se expone el análisis de las descargas del acuífero Ciego- 
Morón, tomando como referencia principal el estudio de Martínez- 
Rodríguez, Hernández-Valdés, Llanusa-Ruiz y  Dilla-Salvador (1986), 
durante el periodo comprendido entre noviembre de 1983 a diciembre de 
1985 (un total de 24 meses). Dichos autores dirigieron su investigación 
hacia el conocimiento del comportamiento hidráulico del mencionado 
acuífero, el cual juega un papel fundamental en el desarrollo 
socioeconómico de la provincia Ciego de Avila, Cuba. El acuífero Ciego- 
Morón se establece en las calizas y dolomías del Mioceno, con espesores 
variables, que pueden alcanzar 300 m, con alta presencia de karst, 
abarcando más de 80% del territorio analizado (Martínez-Rodríguez et al., 
1986). 


Entre los paquetes de programas para realizar el estudio del 
comportamiento de este acuífero se reporta el modelo MEDA (Iglesias, 
1984; Arenillas et al., s.f.), con el objetivo de completar la serie histórica 
de descarga del acuífero en el periodo de calibración seleccionado (de 
noviembre de 1983 a diciembre de 1985). Al no darse detalles del proceso 
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de calibración por parte de Martínez-Rodríguez et a/. (1986) es de suponer 
que ésta fue realizada de forma manual (mediante ensayo y error), pues no 
existe evidencia que indique si la versión utilizada por estos autores es 
aquella que hace uso del algoritmo de Levenberg-Marquardt para la 
calibración automática de sus parámetros, implementada por Erloza 


(1985), citado por Arenillas et al. (s.f.). 


Definición del área a modelar y datos para la estimación de la 
descarga 


En el área en estudio (ver Figura 4), la descarga del acuífero Ciego-Morón 
se produce hacia el norte. Los caudales medidos se reportan en la estación 
de aforo de la cuenca superficial La Yana (área de 1 503 km?) en el periodo 


comprendido en la sección de referencia (Figura 5b). 
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Figura 5. a) Registro de precipitaciones pluviómetro 124. Precipitaciones 
entre diciembre de 1983 y diciembre de 1985, intervalo estudiado por 
Martínez Rodríguez et al. (1986); b) descargas en la estación de aforo de 
La Yana, de diciembre de 1983 a diciembre de 1985. 


Se seleccionó el periodo de precipitaciones entre enero de 1929 a 
diciembre de 2010 (Figura 5a), a partir de los datos de un pluviómetro 
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ubicado en la zona central del acuífero e identificado como P-124, 
coincidiendo con el utilizado por Martínez-Rodríguez et al. (1986), pero 
para una serie histórica entre enero de 1967 y diciembre de 1986. De 
acuerdo con dicho autor, su selección se justificó al resultar imposible hallar 
la isoyeta media mensual, la cual hubiera sido de mucha utilidad para 
establecer la correlación entre la precipitación mensual y los caudales de 
descarga observados en la estación de aforo La Yana. 


Es válido destacar que en dicho estudio se introdujo el volumen de 
descarga de acuerdo con una relación del tipo ecuación (1a); además, no 
se reportaron datos de bombeo y la precipitación útil se toma como el 
acumulado mensual, es decir, sin considerar la evapotranspiración real, 
estableciendo de esta manera los datos de entrada en los modelos 
abordados en esta investigación. 


Experimento numérico 


Un experimento numérico fue conducido con cada uno de los modelos en 
estudio, donde la dimensión del problema quedó previamente establecida a 
partir del número de parámetros, al ser estimados en cada uno de ellos y 
sus posibles variantes. En el algoritmo CS-FMIN, la tolerancia se estableció 
por etapas. Para determinarla se condujeron varios experimentos previos, 
tomando en cuenta que se está empleando una versión débil de CS para 
generar la aproximación inicial a utilizar en FMIN. La tolerancia en CS está 
definida como el mínimo valor que puede obtener la función objetivo en 
estudio, dada en este caso por la ecuación (21), tal y como se muestra en 
la Tabla 3. 


Tabla 3. Experimento numérico. 


Modelo Tolerancia CS Na Pa 

Ecuación (7) 

> 3.2* 5 0.075 
Ecuación (10) 

7 16* 5 0.075 
Ecuación (14) 

5,5** 25 0.25 

% 3.2* 5 0.075 
Ecuación (20) 

5,5** 25 0.25 


(*) Considerando parámetros constantes. 


(**) Considerando variabilidad temporal del parámetro de agotamiento (p. ej., ecuación 
(14): f y ecuación (20): «. 
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En las ejecuciones preliminares de CS-FMIN se detectó que cuando se 
emplea el modelo basado en la curva de agotamiento propuesta por 
Forkasiewicz y Paloc, según la ecuación (14) y considerando el parámetro f 
como constante, la tolerancia de CS debe ser elevada para arribar a una 
aproximación inicial que permita la ejecución de FMIN. Esto no sucede en el 
caso del modelo de descarga no lineal propuesto por el autor, según la 
ecuación (20), donde la tolerancia se mantuvo en los mismos valores de los 
modelos dados por las ecuaciones (7) y (10), lo cual es favorable para esta 
formulación. 


Cuando los parámetros B y xk son considerados variables, tomando en 
cuenta que la dimensión del problema de estudio crece al incrementarse el 
número de parámetros a identificar, el autor prefirió ser conservador y 
admitir una tolerancia ligeramente superior (Tol = 5.5), para agilizar el 
proceso de búsqueda de la aproximación inicial. Por el contrario, la 
tolerancia en la verificación de las restricciones en FMIN se mantuvo 
constante e igual a 1e-06 para todos los casos analizados. 


De igual forma, el número de nidos disponibles Ny y la probabilidad p, de 
que el huevo depositado por el cuco fuera descubierto por otras especies de 
aves se seleccionó a partir de las ejecuciones preliminares a las que se 
hacía referencia con anterioridad. En efecto, pudo detectarse que resulta 
ventajoso que N¿ sea del orden del número de parámetros a ser estimados; 
por esta razón se escogieron los valores de 5 y 25. Cuando el número de 
parámetros se elevó de manera considerable, como consecuencia de 
asignar la variabilidad temporal a f y x«, también fue necesario incrementar 
la probabilidad pa, resultando los valores de 0.075 y 0.25, respectivamente 
(ver Tabla 3). 


Se estableció un total de 20 réplicas para cada modelo en estudio, 
resultando un total de 120 ejecuciones de CS-FMIN, las cuales fueron 
llevadas a cabo en una computadora con 12 Gb de memoria RAM y 
procesador Intel (R) Core (TM) ¡7-4500U CPUO 2.40 GHz. En la Tabla 4 se 
muestra el dominio que representa el espacio de búsqueda para cada uno 
de los parámetros a identificar en cada uno de los modelos, los cuales 
también resultaron de las pruebas preliminares ya mencionadas. 


Tabla 4. Experimento numérico para los modelos en estudio. 


Modelo Dimensión Dominio 
Ecuación (7) m = (1, 70) 
Ecuación (10) d A) 
a = (0, 1e-06) 
m = (1,70) 
Ecuación (14) 3 n=(1, 3) 


p = (0, 2.3e-05)* 


76 
Tecnología y ciencias del agua, 9(5), 58-84, DOI:10.24850/j-tyca-2018-05-03 


Tecnología y e 


CienciaszAgua 


m = (1,70) 
n=(1,3) 
BP; = LBminr BmaxF ¡7 1= 1,...,24 ** 
Praía.= 0, Pro = 2.36-05 
m= (1,70) 
n=(1,3) 
xk = (0, 1e-06)* 
y = (0, 1.0) 
Ecuación (20) m = (1,70) 
n=(1,3) 
27 Ki = X[Kminr Kmáxr ¡7 = 1,...,24 ** 
Kmín = 0, Kmáx = 1e-06 
n = (0, 1.0) 
0; (4%): Ídem Tabla 3, en este particular 24 representa el número de meses. 


26 


Resultados y discusión 


En la Tabla 5 se muestran los parámetros que, en las réplicas conducidas 
por el autor, presentaron el mejor desempeño en la evaluación de la 
función objetivo dada por la ecuación (21). En la Figura 6 se expone el 
comportamiento de los parámetros KB y x, donde puede apreciarse que para 
este último (x) existe una alta correlación entre su valor medio mensual y 
la distribución que minimiza la función objetivo en su mejor ejecución. En 
la Figura 7 se muestra la comparación de caudales de descarga simulados y 
observados para los mejores desempeños de cada modelo. En la Tabla 6 se 
presenta el procesamiento estadístico (media ue, desviación estándar os y 
coeficiente de variación Ce) de las 20 réplicas ejecutadas, donde se detecta 
que el parámetro m exhibe las mayores variaciones con relación a su valor 
medio. 


Tabla 5. Parámetros con mejor desempeño en la evaluación de la F.O. 


Modelo m n a B (e- K n F.O 
- 06) L (ecuación 

(e-07) (e-07) 21) 
Iglesias 68.908 | 2.137 2.86 - a - 3.051 
(1984) 
Tisson et al. 63.341 | 2.247 1.74 - - - 2.871 
(Pérez- 
Franco, 
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1995) 


Forkasiewicz 79.919 | 1.787 - 5.671 - - 15.962 
y Paloc * 
(Pérez- 
Franco, 
1995) 


Forkasiewicz 76.050 1.975 - Ver - - 1.148 
y Paloc ** Figura 


(Pérez- q 


Franco, 
1995) 


Autor * 74.814 | 2.126 - - 2.626 0.809 2.829 


Autor** 67.208 | 2.144 - y Ver 0.788 0.989 
Figura 


7b 


a) b) 


Figura 6. Variabilidad temporal de los parámetros de decrecimiento/agotamiento: a) 
Modelo de Forkasiewicz y Paloc; b) modelo propuesto por el autor. 
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Tabla 6. Resultados estadísticos para los experimentos numéricos. 
Modelos de descarga de acuíferos 
Tisson et al Forkasiewicz- 
Paráme Iglesias ] Paloc* Forkasiewicz-Paloc ** 
(Pérez-Franco á á Autor * Autor ** 
tros (1984) y (Pérez-Franco, (Pérez-Franco, 1995) 
1995) 
1995) 
Mo 0 Cue Mo 0 Cvo Mo O0 Cve Mo O0 Cv Mo O0 Cv Ho 0 Cve 
48.6 | 242. | 0.3 | 42.7 | 283. | 0.3 | 68.0 | 19.2 | 0.0 101.22 57.1 | 184. | 0.2 253.96 
dl 38 1111351 20 36 [959 ¡94 | 20] 75 651%. 8 EN E A a is 
1.87 
2.21 | 0.00 | 0.0 | 2.33 0.00 | 0.0 | 1.81 9 0.0 2.18 | 0.00 | 0.0 
n > 4 29 > 7 37 7 08 1.975 0.001 0.020 0 > 21 2.200 0.006 0.034 
E-04 
2.85 | 4.88 1.72 | 4.05 
ó 3 4 0.0 7 e 0.0 ' - - - - - - - - - - - 
08 12 
E-07 | E-18 E-07 | E-18 
5750 1:01 ver 
B Ñ - y - - - 1 9 0.0 ver ver figura - - E ' y - 
56 | figura 7 | figura 9 
E-06 | E-13 A a 10 
2.57 | 6.16 ver 
a - - - - - - - - - - - - 6 6 0.0 ver ver figura 
30 | figura 7 | figura 9 
E-07 | E-17 z Ñ 10 
0.80 | 0.00 | 0.0 6.707 1.041 
- - - - - - - - - - - - 0.787 
él 3 1 | 37 E-05 E-02 
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Figura 7. Comparación de los caudales simulados y observados: a) Iglesias (1984); b) Tisson et al. (Pérez-Franco, 1995); c) Forkasiewicz y Paloc 

(Pérez-Franco, 1995), con parámetro de decrecimiento f constante; d) Forkasiewicz y Paloc (Pérez-Franco, 1995), con parámetro de decrecimiento 

B variable; e) modelo no lineal propuesto por el autor, con parámetro de agotamiento « constante, f) modelo no lineal propuesto por el autor, con 
parámetro de agotamiento « variable. 
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Tal y como se muestra en la Figura 8, el modelo basado en la ecuación 
de agotamiento dada por Forkasiewicz y Paloc presenta variaciones 
notables del parámetro de decrecimiento de caudal f, lo que redunda en 
una pobre estimación de éste en el experimento conducido, a pesar de 
haber mostrado el mejor desempeño en una de las réplicas de las 20 
ejecutadas. El modelo no lineal presentado por el autor fue muy estable 
en cualquiera de sus alternativas (x constante o variable) y competitivo, 
en comparación con el modelo de Iglesias (1984) y Tisson et al. (Pérez- 
Franco, 1995). 
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Figura 8. Varianza de los parámetros de decrecimiento/agotamiento: a) 
en el modelo de Forkasiewicz y Paloc (Pérez-Franco, 1995); b) en el 
modelo propuesto. 


A pesar de disponer de una sola estación pluviométrica en este estudio 
para validar cada uno de los modelos tratados, es posible llevar a cabo 
simulaciones con cada uno de ellos, con el objetivo de restituir los datos 


1 
Tecnología y ciencias del agua, 9(5), 163-191, DOI:10.24850/j-tyca-2018-03 


Tecnología y 


CienciaszAgua 


históricos dentro del periodo enero de 1929 a diciembre de 2010 (serie 
de 82 años), siendo válido destacar que el modelo de Forkasiewicz y 
Paloc (Pérez-Franco, 1995) se ha excluido en cualquiera de sus 
variantes. En la Figura 8 puede apreciarse una adecuada convergencia 
entre ellos para el régimen de recarga por precipitación a que fue 
sometido el acuífero Ciego-Morón en el periodo de referencia. Cabe 
señalar que aunque se observan descargas máximas mensuales que 
pueden llegar a superar los 10 m*/s, la tendencia de éstas es 
decreciente. 


Conclusiones 


En el presente artículo diferentes modelos agregados para estimar las 
descargas de acuíferos se implementaron en el asistente matemático 
MatLab y acoplados al algoritmo CS-FMIN, demostrándose la robustez 
de la estrategia abordada para reproducir series de caudales de 
descarga en acuíferos kársticos. En cada uno de los modelos se 
identificaron los parámetros relacionados con el volumen de 
precipitación drenado por el acuífero y aquellos que caracterizan el 
agotamiento. 


Fue desarrollado por el autor un modelo de descarga no lineal y su 
desempeño demostró ser satisfactorio en los experimentos numéricos 
conducidos, al ser contrastado con las observaciones en la estación La 
Yana, y con las formulaciones de Iglesias (1984) y Tisson et al. (Pérez- 
Franco, 1995), así como de Forkasiewicz y Paloc (Pérez-Franco, 1995), 
siendo este último desarrollado por sus autores para medios kársticos, 
como el caso aquí abordado. Sin embargo, la demostrada variabilidad 
temporal del parámetro de decrecimiento del caudal limita esta 
formulación sólo al completamiento de series de caudales de descarga 
dentro de intervalos de tiempo específicos y no restituir datos fuera de 
éstos. En este sentido, la propuesta no lineal mostró gran adaptabilidad, 
al considerar parámetros constantes o variables en el tiempo, así como 
las correlaciones a las series observadas, que fueron elevadas. 
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